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Abstract 

A case study in bifurcation and stability analysis is presented, in which reduced dy- 
namical system modelling yields substantial new global and predictive information about 
the behaviour of a complex system. The first smooth pathway, free of pathological and 
persistent degenerate singularities, is surveyed through the parameter space of a nonlin- 
ear dynamical model for a gradient-driven, turbulence-shear flow energetics in magnetized 
fusion plasmas. Along the route various obstacles and features are identified and treated 
appropriately. An organizing centre of low codimension is shown to be robust, several 
trapped singularities are found and released, and domains of hysteresis, threefold stable 
equilibria, and limit cycles are mapped. Characterization of this rich dynamical landscape 
achieves unification of previous disparate models for plasma confinement transitions, sup- 
plies valuable intelligence on the big issue of shear flow suppression of turbulence, and 
suggests targeted experimental design, control and optimization strategies. 

1. Introduction 

The collective dynamics of complex distributed systems often can be usefully described in terms 
of a superposition of rate processes or frequencies which determine the changes in macroscop- 
ically measurable variables as energy flows through the system; that is, a dynamical model 
expressed as a system of coupled ordinary differential equations in a few averaged state vari- 
ables or mode coefficients and several, independently tunable, parameters that represent phys- 
ical properties or external controls. This type of reduced (or low-order or low-dimensional) 
modelling averages over space, mode spectrum structure, single-particle dynamics and other 
details, but the payoff lies in its amenity to sophisticated analytic theory and methods that en- 
able us to track important qualitative features in the collective dynamics, such as singularities, 
bifurcations, and stability changes, broadly over the parameter space. 

Motivated by the need for improved guidance and control of the (mostly bad) behaviour 
of fusion plasmas in magnetic containers, I elaborate in this work a case study in bifurcation 
and stability analysis in which reduced dynamical system modelling yields new global and 
predictive information about gradient driven turbulence-flow energetics that is complementary 
to direct numerical simulation and can guide experimental design. 



Complex 2004 



The 7th Asia-Pacific Conference on Complex Systems 



Reduced dynamical models are powerful tools for describing and analysing complex sys- 
tems such as turbulent plasmas and fluids, primarily because they are supported by well-developed 
mathematics that gives qualitative and global insight, such as singularity, bifurcation, stability, 
and symmetry theory. In principle one can map analytically the bifurcation structure of the en- 
tire phase and parameter space of a reduced dynamical system, but this feat is not possible for 
an infinite-dimensional system, or partial differential equations, and not practicable for systems 
of high order 1 . The usefulness of such models seems to be no coincidence, too: in turbulent 
systems generally, which in detail are both complex and complicated, the dynamics seems to 
take place in a low-dimensional subspace (Holmes et al., 1996). 

It seems paradoxical that enthusiasm for low-dimensional modelling and qualitative analysis 
of fluid and plasma systems has paced the ever larger direct numerical simulations of their flow 
fields. This is an exemplar of how the simplexity and complicity (Cohen and Stewart, 1994) 
juxtaposition can work well: these methods affirm each other, for both have common ground 
in the universal conservation equations for fluid flow (as well as separate bases in mathematics 
and computational science). Developments in one feed developments in the other. Reduced 
dynamical models can give insights into the physics and dynamics of a system in a way that is 
complementary to brute-force numerical simulations of the detailed, spatially distributed mod- 
els from which they are derived. In practice this complementarity means that low-order models 
(which capture few or no spatial modes) can be used to channel information gleaned from the 
generic, qualitative structure of the parameter space — attractors, critical points of onset, sta- 
bility properties, and so on — to numerical simulations (which contain all spatial modes but, on 
their own, bring little physical understanding), giving them purpose and meaning. In turn the 
fluid simulations serve as virtual experiments to validate the low-order approach. 

It is reasonable, therefore, to assert that improved low-dimensional dynamical models for 
plasmas and fluids could provide numerical experimenters with new and interesting challenges 
that will continue to push the limits of computational science and technology. 

/./ Two-dimensional fluid motion in plasmas 

Fusion plasmas in magnetic containers, such as those in tokamak or stellarator experiments, 
are strongly driven nonequilibrium systems in which the kinetic energy of small-scale turbu- 
lent fluctuations can drive the formation of large-scale coherent structures such as shear and 
zonal flows. This inherent tendency to self-organize is a striking characteristic of flows where 
Lagrangian fluid elements see a predominantly two-dimensional velocity field, and is a con- 
sequence of the inverse energy cascade (Kraichnan and Montgomery, 1980). The distinctive 
properties of quasi two-dimensional fluid motion are the basis of natural phenomena such as 
zonal and coherent structuring of planetary flows, but are generally under-exploited in technol- 
ogy- 

In plasmas the most potentially useful effect of two-dimensional fluid motion is suppression 
of high wavenumber turbulence that generates cross-field transport fluxes and degrades con- 
finement (Terry, 2000). Suppression of turbulent transport can manifest temporally as a sponta- 
neous and more-or-less abrupt enhancement of sheared poloidal or zonal flows and concomitant 

'By "reduced" I mean having a phase space dimension of less than four or five, the most influential realm of 
qualitative theories. 
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damping of density fluctuations, and spatially as the rapid development of a localized transport 
barrier or steep density gradient. The phenomenon is often called low- to high-confinement 
(L-H) transitions and has been the subject of intensive experimental, in numero, and theoretical 
and modelling investigations since the 1980s. 

1.2 Confinement transitions and reduced dynamical systems 

The large and lively primary literature on reduced dynamical models for confinement transitions 
and associated oscillations in plasmas represents a sort of consensus on the philosophy behind 
qualitative analysis, if not on the details of the models themselves. What motivates this approach 
is the predictive power that a unified, low-order description of the macroscopic dynamics would 
have in the management of confinement states. Since it is widely acknowledged that control of 
turbulent transport is crucial to the success of the world-wide fusion energy program (U.S. 
DoE, 2000; Fujisawa et al., 2000) it is important to develop predictive models for efficient 
management of access to, and sustainment of, high confinement regimes. For example, if one 
plans to maintain a high confinement state at a relatively low power input by exploiting the 
hysteresis in the transition it would be useful, not to mention cheaper, to know in advance 
which parameters control the shape and extent of hysteresis, or whether it can exist at all in the 
operating space of a particular system, or whether a transition will be oscillatory. 

However, it has been shown (Ball and Dewar, 2000; Ball and Dewar, 2001; Ball et al., 
2002) that many of the models in the literature are structurally flawed. They often contain 
pathological or persistent degenerate (higher order) singularities. An associated issue is that of 
overdetermination, where near a persistent degenerate singularity there may be more defining 
equations than variables. Consequently much of the discussion in the literature concerning 
confinement transitions is qualitatively wrong. Such models cannot possibly have predictive 
power. 

The heart of the matter lies in the mapping between the bifurcation structure and stability 
properties of a dynamical model and the physics of the process it is supposed to represent: if 
we probe this relationship we find that degenerate singularities ought to correspond to some 
essential physics (such as fulfilling a symmetry -breaking imperative, or the onset of hysteresis), 
or they are pathological. In the first case we can usually unfold the singularity in a physically 
meaningful way; in the other case we know that something is amiss and we should revise our 
assumptions. Degenerate singularities are good because they provide opportunities to improve 
a model and its predictive capabilities, but bad when they are not recognized as such. 



The literature on confinement transitions has two basic strands: (1) Transitions are an internal, 
quasi two-dimensional flow, phenomenon and occur spontaneously when the rate of upscale 
transfer of kinetic energy from turbulence to shear and zonal flows exceeds the nonlinear dis- 
sipation rate (Diamond et al., 1992); (2) Transitions are due to nonambipolar ion orbit losses 
near the plasma edge, the resulting electric field providing a torque which drives the poloidal 
shear flow nonlinearly (Itoh and Itoh, 1988; Shaing and Crume, 1989; Stringer, 1993). These 
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two different views of the physics behind confinement transitions are smoothly reconciled for 
the first time in this work. 

A systematic methodology for characterizing the equilibria of dynamical systems involves 
finding and classifying high-order singularities then perturbing around them to explore and map 
the bifurcation landscape (Golubitsky and Schaeffer, 1985). Broadly, this paper is about apply- 
ing singularity theory as a diagnostic tool while an impasto picture of confinement transition 
dynamics is compounded. The bare-bones model is presented in section 2. In section 3. the 
global consequences of local symmetry-breaking are explored, leading to the discovery of an 
organizing centre and trapped degenerate singularities. This leads in to section 4. where I unfold 
a trapped singularity smoothly by introducing another layer that models the neglected physics 
of downscale energy transfer. Section 5. follows the qualitative changes to the bifurcation and 
stability structure that are due to potential energy dissipative losses. In section 6. the unified 
model is presented, in which is included a direct channel between gradient potential energy and 
shear flow kinetic energy. The results and conclusions are summarized in section 7. 

2. The basic model comprises three energy subsystems 

In the edge region of a plasma confinement experiment such as a tokamak or stellarator poten- 
tial energy is stored in a steep pressure gradient which is fed by a power source near the centre. 
Gradient potential energy P is converted to turbulent kinetic energy N, which is drawn off into 
stable shear flows, with kinetic energy F, and dissipation channels. The energetics of this sim- 
plest picture of confinement transition dynamics are schematized in Fig. 1(a). (Nomenclature 
for the quantities here and in the rest of the paper is defined in Table 1 .) 

A skeleton dynamical system for this overall process can be written down directly from 
Fig. 1(a) by inspection: 



The power input Q is assumed constant and the energy transfer and dissipation rates generally 
may be functions of the energy variables. A more physics-based derivation of this system was 
outlined in Ball (2002), in which averaged energy integrals were taken of momentum and pres- 
sure convection equations in slab geometry, using an electrostatic approximation to eliminate 
the dynamics of the magnetic field energy. Equations 1 are fleshed out by substituting specific 



dP_ 

~dt 
dN 

~dt 
dF 

~dt 



1 (P,N,F)-a(P,N,F)-f3(P,N,F) 



Q- 7 {P,N,F) 



a(P,N,F)-n(P,N,F). 



(1) 
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Figure 1 . Energy transfer diagrams for the gradient-driven plasma turbulence-shear flow 
system. Annotated arrows denote rate processes; curly arrows indicate dissipative channels, 
straight arrows indicate inputs and transfer channels between the energy-containing subsys- 
tems. See text for explanations of each subfigure. 



rate-laws for the general rate expressions on the right hand sides: 

dP 
£—— 

dt 

= 1 NP - av ,2 N - (3N 2 

dt 



Q-jNP 



dv' 
i 

' dt 



av'N - fi(P, N)v' 



fi(P, N) =bP- 3/2 + aPN, 



(2a) 

(2b) 

(2c) 
(2d) 



where F = ±v' 2 . The rate expressions in Eqs 2 were derived in Sugama and Horton (1995) and 
Ball (2002) from semi-empirical arguments or given as ansatzes. (Rate-laws for bulk dynamical 
processes are not usually derivable purely from theory, and ultimately must be tested against 
experimental evidence.) 

The rest of this paper is concerned with the character of the equilibria of Eqs 2 and mod- 
ifications and extensions to this system. We shall study the type, multiplicity, and stability of 
attractors, interrogate degenerate or pathological singularities where they appear, and classify 
and map the bifurcation structure of the system. In doing this we shall attempt to answer ques- 
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tions such as: Are Eqs 2 or modified versions a good — that is, predictive — model of the 
system? Does the model adequately reflect the known phenomenology of confinement tran- 
sitions in fusion plasmas? What is the relationship between the bifurcation properties of the 
model and the physics of confinement transitions? 

3. Symmetry-breaking has local and global consequences 

3.1 How to dissolve a pitchfork 



_ip__ 




1 10 100 0.4 

Q 



(b) 



0.6 0.8 

Q 




1 o 1 



Figure 2. (a) (3 = 1, (b) (3 = 50, (b) (3 = (3 crit « 18.58. 7 = I, b = I, a = 0.3, a = 2.4, 
e = 1.5. 



The equilibrium solutions of Eqs 2 are shown in the bifurcation diagrams of Fig. 2, where 
the shear flow v' is chosen as the state variable and the power input Q is chosen as the principal 
bifurcation or control parameter. (In these and subsequent bifurcation diagrams stable equilibria 
are indicated by solid lines and unstable equilibria are indicated by dashed lines.) Several 
bifurcation or singular points are evident. 

The four points in (a) annotated by asterisks, where the stability of solutions changes, are 
Hopf bifurcations to limit cycles, which are discussed in section 3.2. 

On the line v' — the singularity P is found to satisfy the defining and non-degeneracy 
conditions for a pitchfork, 

G = G c = G (( = G x = 0, G a( + 0, G C x + 0, (3) 

where G is the bifurcation equation derived from the zeros of Eqs 2, ( represents the chosen state 
variable, A represents the chosen control or principal bifurcation parameter, and the subscripts 
denote partial derivatives. In the qualitatively different bifurcation diagrams (a) and (b) the 
dissipative parameter f3 is relaxed either side of the critical value given in (c), where the perfect, 
twice-degenerate pitchfork is represented. Thus for (3 < (3 CI i t (a poorly dissipative system) the 
turning points in (a) appear and the system may also show oscillatory behaviour. The dynamics 
are less interesting for [3 > f3 crit (a highly dissipative system) as in (b) because the turning 
points, and perhaps also the Hopf bifurcations, cannot occur. 

However, P is persistent through variations in [3 or any other parameter in Eqs 2. (This fact 
was not recognized in some previous models for confinement transitions, where such points 
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were wrongly claimed to represent second-order phase transitions.) Typically the pitchfork 
is associated with a fragile symmetry in the dynamics of the modelled physical system. The 
symmetry in this case is obvious from Fig. 2: in principle the shear flow can be in either 
direction equally. In real life (or in numero), experiments are always subject to perturbations 
that determine a preferred direction for the shear flow, and the pitchfork is inevitably dissolved. 
In this case the perturbation is an effective force or torque from any asymmetric shear-inducing 
mechanism, such as friction with neutrals in the plasma or external sources, and acts as a shear 
flow driving rate. Assuming this rate to be small and independent of the variables over the 
characteristic timescales for the other rate processes in the system, we may revise the shear 
flow evolution Eq. 2c as 

dv' 

2— = av'N - fi(P, N)v' + <p, (4) 

where the symmetry-breaking term (p models the shear flow drive. The corresponding energy 
transfer schematic is shown in Fig. 1(b). 

The pitchfork P in Fig. 2 (c) can now be obtained exactly by applying the conditions (3) to 
the zeros of Eqs 2a, 2b, and 4, with 2d, and with ( = v' and A = Q: 

(ct 2 i 2 2a 3, j\/a/a \ 

The other singularity T on v ' = satisfies the defining and non-degeneracy conditions for a 
transcritical bifurcation, 

G = G c = G x = 0, G cc ^ 0, det [ ^ cc ^ AC ) = det d 2 G < 0. (5) 

\(jac G\\J 

It is once-degenerate and also requires the symmetry-breaking parameter for exact definition. 



3.2 A walk along untrodden ways 

A bifurcation diagram where P is fully unfolded, that is, for ip ^ and f3 ^ f3 crit , is shown in 
Fig. 3. This diagram is rich with information that speaks of the known and predicted dynamics 
of the system and of ways in which the model can be improved further, and which cannot be 
inferred or detected from the degenerate bifurcation diagrams of Fig. 2. It is worthwhile to step 
through Fig. 3 in detail, with the energy schema Fig. 1 (b) at hand. Let us begin on the stable 
branch at Q ~ 0.1. 

Here the pressure gradient is being charged up, the gradient potential energy is feeding the 
turbulence, and the shear flow is small but positive because the sign of the perturbation if is 
positive. As the power input Q is increased quasistatically the shear flow begins to grow but 
at the turning point, where solutions become unstable, there is a discontinuous transition to the 
upper stable branch in v' and P, (a) and (b), and to the lower stable branch in N, (c). At the 
given value of f3 the hysteresis is evident: if we backtrack a little the back transition takes place 
at a lower value of Q. Continuing along the upper stable branch in v' we encounter another 
switch in stability; this time at a Hopf bifurcation to stable period one limit cycles. (In this and 
other diagrams the amplitude envelopes of limit cycle branches are marked by large solid dots.) 
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Figure 3. = 0.05, other parameters as for Fig. 2 (a). 



From the amplitude envelope we see that the oscillations grow as power is fed to the system 
then are extinguished rather abruptly, at a second Hopf bifurcation where the solutions regain 
stability. The shear flow decreases toward zero as the pressure-dependent anomalous viscosity, 
the second term in Eq. 2d, takes over to dissipate the energy at high power input. 

The system may also be evolved to an equilibrium on the antisymmetric, — v' branch in 
Fig. 3 (a), by choosing initial conditions 2 appropriately or a large enough kick. However, if the 
power input then falls below the turning point at (V, Q) xs (—0.9, 0.46) we see an interesting 
phenomenon: the shear flow spontaneously reverses direction. The transient would nominally 
take the system toward the nearest stable attractor, the lower +v' branch, but since it would then 
be sitting very close to the lower +v' turning point small stochastic fluctuations could easily 
induce the transition to the higher +v' branch. See the inset zoom-in over this region in (a). 
Here is an example of a feature that is unusual in bifurcation landscapes, a domain over which 
there is fivefold multiplicity comprising three stable and two unstable equilibria. Two more 
examples of threefold stable domains will be shown in section 5. 

The same equilibria depicted using P and N as dynamical variables in (b) and (c) are anno- 
tated to indicate whether they correspond to the +v' or —v' domain. For clarity the amplitude 
envelopes of the limit cycle solutions are omitted from (b) and (c). 

In the remainder of this paper I concentrate on the +v' branches and ignore the —v' domain. 

2 A note on initial conditions: The stable attractor that is spontaneously sought by a dynamical system will be 
the one for which the initial conditions lie within its basin of attraction. The overall bifurcation structure of the 
type of dynamical system considered here is independent of initial conditions. 



Complex 2004 



The 7th Asia-Pacific Conference on Complex Systems 



3.3 A metamorphosis of the dynamics 

Now we approach the very heart of the model, the organizing centre; strangely enough via the 
branch of unstable solutions that is just evident in Fig. 3 (a) and (b) in the top left-hand corner 
and (c) in the lower left hand corner. The effects of symmetry-breaking are more far-reaching 
than merely providing a local universal unfolding of the pitchfork, for this branch of equilibria 
was trapped as a singularity at (V, Q) = (oo, 0) for uo = 0. The organizing centre itself, 
described as a metamorphosis in Ball (2002), can be encountered by varying (p. The sequence 
in Fig. 4 tells the story visually. 




0.1 1 10 100 0.1 1 10 100 

Q Q 

Figure 4. (a) ip = 0.08, (b) up « 0.08059 = <p Tm , (c) up = 0.1, (d) up = 0.11. e = 1, other 
parameters as for Fig. 2 (a). 

The "new" unstable branch develops a Hopf bifurcation at the turning point. (Strictly speak- 
ing, this is a degenerate Hopf bifurcation, called DZE, where a pair of complex conjugate eigen- 
values have zero real and imaginary components.) As uo is tuned up to 0.08 (a) a segment of 
stable solutions becomes apparent as the Hopf bifurcation moves away from the turning point; 
the associated small branch of limit cycles can also just be seen. At ud T m, the metamorphosis, 
the "new" and "old" branches exchange arms, (b). 

The metamorphosis satisfies the conditions (5) and is therefore an unusual, non- symmetric, 
transcritical bifurcation. It signals a profound change in the type of dynamics that the system 
is capable of. For ip > uo TM a transition must still occur at the lower limit point, but there is 
no classical hysteresis, (c) and (d). In fact classical hysteresis is (locally) forbidden by the non- 
degeneracy condition ^ in Eqs 5. Various scenarios are possible in this regime, including 



Complex 2004 



The 7th Asia-Pacific Conference on Complex Systems 



a completely non-hysteretic transition, a forward transition to a stable steady state and a back 
transition from a large period limit cycle, or forward and back transitions occurring to and from 
a limit cycle. 

3.4 The story so far 

The symmetry-broken model comprises Eqs 2a, 2b, and 4, with 2d. The bifurcation structure, 
some of which is depicted in Figs 3 and 4, predicts various behaviours: 

• Shear flow suppression of turbulence; 

• Smooth, hysteretic, non-hysteretic, and oscillatory transitions; 

• Spontaneous and kicked reversals in direction of shear flow; 

• Saturation then decrease of the shear flow with power input due to pressure-dependent 
anomalous viscosity; 

• A metamorphosis of the dynamics through a transcritical bifurcation. 

A critical appraisal of experimental evidence that supports the qualitative structure of this model 
is given in Ball (2004). With the exception of the last item all of the above dynamics have been 
observed in magnetically contained fusion plasma systems. The model would therefore seem 
to be a "good" and "complete" one, in the sense of being free of pathological or persistent 
degenerate singularities and reflecting observed behaviours. 

However, there are several outstanding issues that suggest the model is still incomplete. One 
issue arises as a gremlin in the bifurcation structure that makes an unphysical prediction, another 
comes from a thermal diffusivity term that was regarded as negligible in previous work on this 
model. A third issue arises from the two strands in literature on the physics of confinement 
transitions: the model as it stands does not describe confinement transitions due to a nonlinear 
electric field drive. 

4. Shear flows also generate turbulence 

The first issue of incompleteness concerns a pathology in the bifurcation structure of the model, 
implying infinite growth of shear flow as the power input falls. Before we pinpoint the culprit 
singularity, it is illuminating to evince the physical — or unphysical — situation through a 
study of the role of the thermal capacitance parameter e, which regulates the contribution of the 
pressure gradient dynamics, Eq. 2a, to the oscillatory dynamics of the system. Conveniently, 
we can use e as a second parameter to examine the stability of steady-state solutions around the 
Hopf bifurcations in Fig. 4 without quantitative change. The machinery for this study consists 
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of the real-time equations 2a, 2b, and 4 recast in "stretched time" r — t/e, 

dP 



e (jNP - av' 2 N - fiN 2 ) 



dN 

~dr~ 

dv' 

2— = e (av'N - /i(P, N)v' + <p) , 
dr 



(6a) 
(6b) 
(6c) 



and the two-parameter locus of Hopf bifurcations in the real-time system shown in Fig. 5. We 
consider two cases. 




Figure 5. The curve is the locus 
of the Hopf bifurcations in Fig. 
4(c) over variations in e. 



1. The high-capacitance regime: 

The maximum on the curve in Fig. 5 marks a degenerate Hopf bifurcation, a DZE point. 
Here the two Hopf bifurcations at lower Q in Fig. 4(c) merge and are extinguished as 
e increases. (This merger through a DZE has obviously occurred in Fig. 4(d), where ip 
is varied rather than e.) The surviving upper Hopf bifurcation moves to higher Q as e 
increases further. 

In this high-capacitance regime the dynamics becomes quasi one-dimensional on the 
stretched timescale. To see this formally, define a — 1/e and multiply the stretched- 
time equations 6b and 6c through by a. In the limit a — > the kinetic energy variables 
are slaved to the pressure gradient (or potential energy) dynamics. Switching back to 
real time and multiplying Eq. 2a through by a, for a — > P « P . The kinetic energy 
subsystems see the potential energy as a constant, "infinite source". It is conjectured that 
the surviving upper Hopf bifurcation moves toward (Q, v') = (oo, oo) and for e ^> 1 the 
dynamics becomes largely oscillatory in real time, with energy simply sloshing back and 
forth between the turbulence and the shear flow. 

2. The low-capacitance regime: 

The minimum on the curve in Fig. 5 marks another degenerate Hopf bifurcation, also 
occurring at a DZE point. Here the two Hopf bifurcations at higher Q in Fig. 4(c) merge 
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and are extinguished as e decreases. The surviving Hopf bifurcation moves to lower Q 
— and higher v' — as e decreases further. This scenario is illustrated in Fig. 6, where 
the steady- state curve and limit cycle envelope are roughly sketched in a decreasing e 
sequence. 



>- decreasing e 

i i i 



\ 

\. 


» 
\ 




v' 


1 


\ 


\ 






\ . 


s '*•••„ 


\ 

\ 



















Q 

Figure 6. As e is decreased further than the minimum in Fig. 5 the remaining Hopf bi- 
furcation slides up the steady-state curve, which becomes stable toward unrealistically 
high v' and low Q. 

In this low-capacitance regime the dynamics also becomes quasi one-dimensional, and as 
e — > the conjectured fate of the surviving Hopf bifurcation is a double zero eigenvalue 
trap at (Q, v') = (0, oo). To see why this can be expected, consider again the stretched- 
time system, Eqs 6. For eClwe have N m N and v' rs v . On the stretched timescale 
the potential energy subsystem sees the kinetic energy subsystems as nearly constant, and 
P Pa (P -Q/ (iV 7)) exp (-7Vo7t) +Q/ (iV 7). Reverting to real time, as edP/dt -> 
we have P m Q/ ('fN); the potential energy is reciprocally slaved to the kinetic energy 
dynamics. 

The anomaly in this low-capacitance picture is that, as the power input Q ebbs, the shear 
flow can grow quite unrealistically. With diminishing e the Hopf bifurcation moves up- 
ward along the curve, the branch of limit cycles shrinks, and the conjugate pair of pure 
imaginary eigenvalues approaches zero. It would seem, therefore, that some important 
physics is still missing from the model. 

4.1 A trapped singularity is found and released 

What is not shown in Figs 3 and 4 (because a log scale is used for illustrative purposes) is a 
highly degenerate branch of equilibria that exists at Q = where N = and v' = (P 3 ^ 2 (p)/b; 
it is shown in Fig. 7(a). For (p > there is a trapped degenerate turning point, annotated as s4, 
where the "new" branch crosses the Q = branch. 

The key to its release (or unfolding) lies in recognizing that kinetic energy in large-scale 
structures inevitably feeds the growth of turbulence at smaller scales, as well as vice versa 
(Terry, 2000). In a flow where Lagrangian fluid elements locally experience a velocity field that 
is predominantly two-dimensional there will be a strong tendency to upscale energy transfer 
(or inverse energy cascade, see Kraichnan and Montgomery (1980)), but the net rate of energy 
transfer to high wavenumber (or Kolmogorov cascade, see Ball (2004)) is not negligible. What 
amounts to an ultraviolet catastrophe in the physics when energy transfer to high wavenumber is 
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Figure 7. (a) k = 0, (b) k = 0.001. up = 0.08, (3 = 0.3, b = 1, a = 0.3, a = 2.4, 7 = 1, e 
77ze starred stability changes are Hopf bifurcations. 



neglected maps to a trapped degenerate singularity in the mathematical structure of the model. 
The trapped singularity s4 may be unfolded smoothly by including a simple, conservative, back- 
transfer rate between the shear flow and turbulent subsystems: 

dN 

— = 1 NP - av' 2 N - (3N 2 + kv' 2 (7a) 
dt 

dv' 

2— = av'N - fi(P, N)v' + ud - kv'. (7b) 

(JjL 

The model now consists of Eqs 7 and 2a, with 2d, and the corresponding energy transfer 
schematic is Fig. 1(c). The back- transfer rate coefficient k need not be identified with any 
particular animal in the zoo of plasma and fluid instabilities, such as the Kelvin-Helmholtz 
instability; it is simply a lumped dimensionless parameter that expresses the inevitability of 
energy transfer to high wavenumber. 

The manner and consequences of release of the turning point s4 can be appreciated from 
Fig. 7(b), from which we learn a salutary lesson: unphysical equilibria and singularities should 
not be ignored. The unfolding of s4 creates a maximum in the shear flow, and (apparently) 
a fourth Hopf bifurcation is released from a trap at infinity. At the given values of the other 
parameters this unfolding of s4 has the effect of forming a finite-area isola of steady-state solu- 
tions, but it is important to visualize this (or, indeed, any other) bifurcation diagram as a slice 
of a three-dimensional surface of steady states, where the third coordinate is another parameter. 
(Isolas of steady- state solutions were first reported in the chemical engineering literature, where 
nonlinear dynamical models typically include a thermal or chemical autocatalytic reaction rate 
(Ball, 1999).) 

In Fig. 8 we see two slices of this surface, prepared in order to demonstrate that the metamor- 
phosis identified in section 3.3 is preserved through the unfolding of s4. Here the other turning 
points are labelled si, s2, and s3. Walking through Fig. 8 we make the forward transition at 
si and progress along this branch through the onset of an a limit cycle regime, as in Fig. 3. 
For obvious reasons we now designate this segment as the intermediate shear flow branch, and 
the isola or peninsula as the high shear flow branch. In (a) a back-transition occurs at s2. The 
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Figure 8. (a) up = 0.083, (b) up = 0.084. k = 0.001, (3 = 0.3, b = 1, a = 0.3, a = 2.4, 7 
e = l. Starred points are Hopf bifurcations. 



1, 



system can only reach a stable attractor on the isola by a transient, either a non-quasistatic jump 
in a second parameter or an evolution from initial conditions within the appropriate basin of 
attraction. In (b) as we make our quasistatic way along the intermediate branch with diminish- 
ing Q the shear flow begins to grow, then passes through a second oscillatory domain before 
reaching a maximum and dropping steeply; the back transition in this case occurs at s4. 



5. Thermal dissipation affects the bifurcation structure 

In the model so far the only outlet channel for the potential energy is conversion to turbulent 
kinetic energy, given by the conservative transfer rate ^PN. However, in a driven dissipative 
system such as a plasma other conduits for gradient potential energy may be significant. 

The cross-field thermal diffusivity, a neoclassical transport quantity (Braginskii, 1965) is 
often assumed to be negligible in the strongly-driven turbulent milieu of a tokamak plasma 
(Sugama and Horton, 1994; Sugama and Horton, 1995; Ball et al., 2002), but here Eq. 2 is 
modified to include explicitly a linear "infinite sink" thermal energy dissipation rate: 

e< ^ = Q- 1 NP ~ XP- (8) 

Following Thyagaraja et al. (1999) x is taken as as a lumped dimensionless parameter and 
the rate term %P as representing all non-turbulent or residual losses such as neoclassical and 
radiative losses. The model now consists of Eqs 7 and 8, with 2d and the corresponding energy 
schematic is Fig. 1(d). 

This simple dissipative term has profound effects on the bifurcation structure of the model, 
and again the best way to appreciate them is through a guided walking tour of the bifurcation 
diagrams. In Fig. 9 the series of bifurcation diagrams has been computed for increasing values 
of x an d a connected slice of the steady state surface (i.e., using a set of values of the other 
parameters for which the metamorphosis has already occurred). 
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Figure 9. (a) X = 0.005, (b) x = 0.01, (c) \ = 0.05, (d) x = 0.1, (e) x = 0.2. « = 0.001, 
</? = 0.088, /3 = 0.3, b = 1, a = 0.3, ct = 2.4, j — 1, e — 1. In (a) to (d) the Hopf bifurcations 
are starred, and for clarity the limit cycle amplitude envelopes are omitted. 



A qualitative change is immediately apparent, which has far-reaching consequences: for 
X > the two new turning points s5 and s6 appear, born from a local cusp singularity that was 
trapped at X = 0. Overall, from (a) to (e) we see that si does not shift significantly but that the 
peninsula becomes more tilted and shifts to higher Q, but let us begin the walk at s 1 in (b). Here, 
as in Fig. 8, the transition occurs to an intermediate shear flow state and further increments of 
Q take the system through an oscillatory regime. But the effect of decreasing Q is radically 
different: at s6 a discontinuous transition occurs to a high shear flow state on the stable segment 
of the peninsula. From this point we may step forward through the shear flow maximum and 
fall back to the intermediate branch at s5. We see that over the range of Q between s5 and s6 
the system has five steady states, comprising three stable interleaved with two unstable steady 
states. As in Fig. 8(b) a back transition at low Q occurs at s4. 

The tristable regime in (b) has disappeared in (c) in a surprisingly mundane way: not through 
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a singularity but merely by a shift of the peninsula toward higher Q. But this shift induces a 
different tristable regime through the creation of s7 and s8 at another local cusp singularity. In 
(d) s4 and s7 have been annihilated at yet another local cusp singularity. It is interesting and 
quite amusing to puzzle over the 2-parameter projection of these turning points si, s8, s7, and 
s4 followed over x, it is given in Fig. 10. The origins of the three local cusps can be read off 
the diagram, keeping in mind that the crossovers are a trompe de I'oeil: they are nonlocal. 




Figure 10. The turning points si, s8, 
s7, and s4 are followed over X- K = 
0.001, up = 0.088, = 0.3, 6=1, 
a = 0.3, a = 2.4, 7 = 1, e = 1. 



At s5 in Fig. 9(c), (d), and (e) the system transits to a limit cycle, rather than to a stable 
intermediate steady state. Shown in Fig. 1 1 are the bifurcation diagrams in N and P corre- 
sponding to Fig. 9(e). The pressure gradient jumps at si because the power input exceeds the 
distribution rates, and oscillatory dynamics between the energy subsystems sets in abruptly at 
s5. The turbulence is enormously suppressed due to uptake of energy by the shear flow, but 
rises again dramatically with this hard onset of oscillations. 




Q Q 
Figure 11. The turbulence and pressure gradient subsystems corresponding to Fig. 9(e). 



6. A unified model includes a nonlinear shear flow drive 



The early theoretical work on confinement transitions attempted to explain edge L-H transitions 
exclusively in terms of the electric field driving torque created by nonambipolar ion orbit losses 
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(Itoh and Itoh, 1988; Shaing and Crume, 1989; Stringer, 1993), with no coupling to the internal 
dynamics of energy transfers from the potential energy reservoir in the pressure gradient. The 
electric field is bistable, hence the transition to a high shear flow, or high confinement, regime 
is discontinuous and hysteretic. Although there are many supporting experiments (Fujisawa, 
2003), this exposition of the physics behind confinement transitions is incomplete because it 
cannot explain shear flow suppression of turbulence — a well-known characteristic of L-H 
transitions. 

Here this "electric field bifurcation" physics is treated as a piece of a more holistic physical 
picture and a simple model for the rate of shear flow generation due to this physics is used 
to create a unified dynamical model for confinement transitions. Following the earlier authors 
this rate is given as r{P) = z/exp — (w 2 /P) 2 , which simply says that the rate at which ions 
are preferentially lost, and hence flow is generated, is proportional to a collision frequency v 
times the fraction of those collisions that result in ions with sufficient energy to escape. The 
form of the energy factor assumes an ion distribution that is approximately Maxwellian and w 2 , 
analogous to an activation energy, is proportional to the square of the critical escape velocity. 
In this form of the rate expression I have explicitly included the temperature-dependence of r, 
through P, which couples it to the rest of the system. If w is high the rate is highly temperature 
(pressure gradient) sensitive. (For heuristic purposes constant density is assumed, constants and 
numerical factors are normalized to 1, and the relatively weak temperature dependence of the 
collision rate v is ignored.) 

For convenience the equations for the unified model are gathered together: 

dP 

e— = Q- 7 iVP - v' 2 r(P) - X P (9a) 

(JjL 

= 7 iVP - av' 2 N - (3N 2 + kv' 2 (7a) 

at 
dv' 

2— = av'N - fi(P, N)v' + v'r(P) - kv' + <p (9b) 

(-lib 



r(P) = z/exp 



w 2 /P" 2 



(9c) 



fi(P, N) = bp- 3/2 + aPN. (2d) 

The corresponding energy schematic is Fig. 1(e) where it is seen that r(P) is a competing 
potential energy conversion channel, that can dominate the dynamics when the critical escape 
velocity w is low or the pressure is high. This is exactly what we see in the bifurcation diagrams, 
Fig. 12. Overall, the effect of this contribution to shear flow generation from the ion orbit 
loss torque is to elongate and flatten the high shear flow peninsula. The Hopf bifurcations 
that are starred in (a), where the contribution is relatively small, have disappeared in (b) at 
a DZE singularity. What this means is that as r{P) begins to take over there is no longer 
a practicably accessible intermediate branch, as can be seen in (c) where the intermediate 
branch is unstable until the remaining Hopf bifurcation is encountered at extremely high Q. 
Locally, in the transition region, as r(P) becomes significant the bifurcation diagram begins 
to look more like the simple S-shaped, cubic normal form schematics with classical hysteresis 
presented by earlier authors. However, this unified model accounts for shear flow suppression 
of the turbulence (d), whereas theirs could not. 
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7. Summary and conclusions 

The generation of stable shear flows in plasmas, and the associated confinement transitions and 
oscillatory behaviour in tokamaks and stellarators, is regulated by Reynolds stress decorrelation 
of gradient-driven turbulence and/or by an induced bistable radial electric field. These two 
mechanisms are smoothly unified by the first smooth road through the singularity and 
bifurcation structure of a reduced dynamical model for this system. 

The model is constructed self-consistently, beginning from simple rate-laws derived from 
the basic pathways for energy transfer from pressure gradient to shear flows. It is iteratively 
strengthened by finding the singularities and allowing them to "speak for themselves", then 
matching up appropriate physics to the unfoldings of the singularities. 

The smooth road from turbulence driven to electric field driven shear flows crosses interest- 
ing territory: 

• Hysteresis is possible in both regimes and is governed by different physics. 

• A metamorphosis of the dynamics is encountered, near which hysteretic transitions are 



Complex 2004 



The 7th Asia-Pacific Conference on Complex Systems 



forbidden. The metamorphosis is a robust organizing centre of codimension 1, even 
though there are singularities of higher codimension in the system. 

• Oscillatory and tristable domains are encountered. 

• To travel the smooth road several obstacles are successively negotiated in physically 
meaningful ways: a pitchfork is dissolved, simultaneously releasing a branch of solu- 
tions from a singular trap at infinity, a singularity is released from a trap at zero power 
input, and another is released from a trap at zero thermal diffusivity. 

In particular, these results suggest strategies for controlling access to high confinement states 
and manipulating oscillatory behaviour in fusion experiments. More generally I have shown that 
low-dimensional models have a useful role to play in the study of one of the most formidable of 
complex systems, a strongly driven turbulent plasma. Having survived such a trial-by-ordeal, 
the methodology is expected to continue to develop as a valuable tool for taming this and other 
complex systems. 



Table 1. Glossary of nomenclature 


P 


pressure gradient potential energy 


N 


kinetic energy of the turbulence 


F 


shear flow kinetic energy 


v' 


shear flow 


e 


dimensionless thermal capacitance 


Q 


power input to the plasma pressure gradient 


7 


turbulence driving rate coefficient 


a 


shear flow driving rate coefficient due to Reynolds stress 




decorrelation of turbulence 


P 


turbulence dissipation rate coefficient 


M 


kinematic viscosity 


b 


neoclassical viscosity 


a 


anomalous viscosity 




symmetry-breaking shear flow perturbation 


K 


energy back-transfer (downscale) rate coefficient 


X 


thermal diffusivity 


V 


ion-ion collision rate 


w 


critical ion escape velocity 
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